home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / SIMPS.LIB < prev    next >
Text File  |  1985-04-03  |  1KB  |  50 lines
  1.  
  2.  
  3. { -> 278 }
  4. procedure simps(function f(x: real): real;
  5.         lower,upper,tol    : real;
  6.         var sum        : real);
  7.  
  8. { numerical integration by Simpson's rule }
  9. { function is f (as parameter), limits are lower and upper }
  10. { with number of regions equal to pieces }
  11. { partition is delta_x, answer is sum }
  12.  
  13. var    i        : integer;
  14.     x,delta_x,even_sum,
  15.     odd_sum,end_sum,
  16.     end_cor,sum1    : real;
  17.     pieces        : integer;
  18.  
  19. function dfx(x:real):real;
  20. begin
  21.   dfx:=-1.0/sqr(x)
  22. end;
  23.  
  24. begin
  25.   pieces:=2;
  26.   delta_x:=(upper-lower)/pieces;
  27.   odd_sum:=f(lower+delta_x);
  28.   even_sum:=0.0;
  29.   end_sum:=f(lower)+f(upper);
  30.   end_cor:=dfx(lower)-dfx(upper);
  31.   sum:=(end_sum+4.0*odd_sum)*delta_x/3.0;
  32.   repeat
  33.     pieces:=pieces*2;
  34.     sum1:=sum;
  35.     delta_x:=(upper-lower)/pieces;
  36.     even_sum:=even_sum+odd_sum;
  37.     odd_sum:=0.0;
  38.     for i:=1 to pieces div 2 do
  39.       begin
  40.     x:=lower+delta_x*(2.0*i-1.0);
  41.     odd_sum:=odd_sum+f(x)
  42.       end;
  43.     sum:=(7.0*end_sum+14.0*even_sum+16.00*odd_sum
  44.                 +end_cor*delta_x)*delta_x/15.0;
  45.   until abs(sum-sum1)<=abs(tol*sum)
  46. end;    { simps }
  47.  
  48.  
  49.  
  50.